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Abstract 

This article deals with the numerical integration in time of nonlinear Schrodinger equations. The main 
application is the numerical simulation of rotating Bose-Einstein condensates. The authors perform a change 
of unknown so that the rotation term disappears and they obtain as a result a nonautonomous nonlinear 
Schrodinger equation. They consider exponential integrators such as exponential Runge-Kutta methods and 
Lawson methods. They provide an analysis of the order of convergence and some preservation properties 
of these methods in a simplified setting and they supplement their results with numerical experiments with 
realistic physical parameters. Moreover, they compare these methods with the classical split-step methods 
applied to the same problem. 
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Keywords. Gross-Pitaevskii equation, nonlinear Schrodinger equation, exponential Runge-Kutta and Lawson methods, 
splitting methods, superconvergence. 


1 Introduction 


This paper deals with the numerical integration of nonautonoumous nonlinear Schrodinger (NLS) equations 
which read 


f idtip =—^Aip + V(t,x)'tp + (t,x) G K+ X 

\ 'tp(0,x) = V’o(x), X e 


The operator A denotes the usual Laplace operator on d G N*. The potential function V is smooth 
and K G N. The unknown is a complex-valued wavefunction associated to the given initial datum bo- The 
derivation and the analysis of efficient semi-discrete numerical methods for the time integration of these equations 
have a long history. Some authors are interested in the finite time accuracy of the schemes [2H1 IMl HI m]. 
With additional hypotheses corresponding to physically relevant situations, equation (1.1) may have several 
invariants such as energy, momentum, etc, ..., in addition to mass conservation. Several authors show interest 
in the preservation of the invariants by the numerical schemes HHUSSISllllS]- Beyond finite time integration 


1 



2 


C. BESSE AND G. DUJARDIN AND 1. LACROIX-VIOLET 


of Eq. 0, asymptotic regimes have been considered: long time preservation of invariants [IIIMIIIH] and 
semi-classical regimes [SI dHl ESI EHl El] • 

Our goal in this paper is to introduce and analyse two classes of exponential methods for the time integration 
of (1.1 1 on a d-dimensional space. The first class is that of exponential Runge-Kutta methods I3Z1IM1E7]. The 
second class relies on the Lawson techniques [131113 . We focus on situations were the dynamics of equation (1.1 1 
essentially stays in a bounded domain of K'^. Hence, we replace Eq. (1.1 1 with the same equation on a large 
periodic torus (with characteristic size <5 > 0). Such a change is usually not uniform with respect to 5 (see 
subsection 1.2 1 . However it is physically relevant even if it is only valid for reasonably bounded times a priori. 
Our main application is the Gross-Pitaevskii equation modeling a rotating Bose-Einstein Condensate (BEG) 
in for which the dynamics of the solutions is much spatialy localized. We prove our results for functions 
in Sobolev spaces for cr > d/2, instead of The fully discrete corresponding situations are 

not similar. Discretizing functions on the whole space requires to deal with carefully chosen specific boundary 
conditions. Working on a torus has the numerical advantage of avoiding boundary conditions. Moreover the 
numerical computation of the linear group generated by zA is much simpler in the case of the torus m- An 
alternative to the reduction from to T/ consists in using pseudospectral methods introduced in |7] and 
analyzed in m- In contrast to the approach presented in this paper, pseudospectral methods can not be used 
for general confining potentials {e.g. quadratic potentials without symmetry in the rotation plane ( 7 ^ 7 ^ in 

(1.3l) that result in non-autonomous linear parts after the change of variable introduced below (see (1.9l) or 
quadratic-quartic potentials such as V{x) = ||a;|p + (1 + sin(||a:||^))||a;||^). 

Before the introduction of the methods (Sections[^and [^ an d our numerical results (Section]^, we introduce 
below our main application and show how it fits equation ( 1 . 11 . 


1.1 Presentation of the application 

A Bose-Einstein Gondensate is the state of matter reached by a dilute gas of bosons cooled to very low temper¬ 
ature. A large fraction of bosons occupies the lowest quantum state so that macroscopic quantum phenomena 
become apparent. This phenomena was theoretically predicted by Bose in 1924 for photons ESI and generalized 
to atoms by Einstein in 1925 jSD]- Tii® first experimental evidence of BEG was obtained in 1995 [T1 123j. 

At low temperature, a rotating planar BEG can be described by the macroscopic complex-valued wave 
function (p = ip(t, x) whose evolution is governed by a Gross-Pitaevskii Equation (GPE) with an angular 
momentum rotation term. After a suitable changes of variables [3j, the dimensionless GPE in d = 2 dimensions 
satisfied by (p can be written for x G 

idtp= -]^/i^p + Vc{yi)p + l3\p\'^p-nRp. ( 1 . 2 ) 

The real-valued function W = Vc{x) corresponds to a smooth potential depending only on the space variables 
denoted by x = {x,yy. In our physical context, this potential is confining: this means that 14(x) tends to -l-oo 
when ||x|| = ^x^ tends to -l-oo. For example, in this paper, we consider potentials of the form 

W(x) = ^ {jlx^ + 'lyy^), (1.3) 


where ^xily > 0- This confining potential competes with the rotation operator —flR = iU,(xdy — yd^) at 
angular speed H G K: the former tends to make bosons stay together at the origin of the plane, while the latter 
tends to spread the bosons out. The real coefficient /3 represents the nonlinearity strength, and comes from the 
averaged effect of the bosons. The evolution equation (1.2) is supplemented with an initial condition 

93(0,x) = y)o(x), for all x G (1.4) 


If one introduces the functional 

E{p) = 


+ Kjpy + ^\py - nRe{p*Rp) 


dx, 


(1.5) 


then ( 1 . 2 ) reads 


idtpit,.) = 

The equation happens to be Hamiltonian, and the energy is preserved by the dynamics: along a solution t 1 —> 


p{t ,.) of (1.2 )-(1.4), one has E{p(t ,.)) = E{po). In addition to preserving the energy, the evolution preserves the 


mass of the wave function: along a solution 1 1 —>■ .) of (1.2)-(1.4), one has /jja |<^(^,x |(/7o(x)pdx. 

Another dynamical feature of this equation is the evolution of the angular momentum expectation: if one 
denotes by 

<R>{t)= / (p*(t,x)i?(^(t,x)dx, (1.6) 

dK2 
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then this real-valued function is constant in the special case of a radial harmonic potential ( 73 ; = 7 ^ in (1.3l) 
and has a more complex dynamics in more general cases (see Lemma 6.2.1 in |48j 1. 

In the last decades this model has been studied a lot [SlEHIlllS]. An important issue in the numerical time 
integration of equations (1.2)-(1.41 comes from the rotation term R. Following [H], we introduce new coordinates 
that allow us to put equation (1.2 1 in the form of equation (l.ll. Let us set for t G K, 


^ _ /cos(fit) — sin(nt)\ 
\^sin(Vlf) cos(nt) j 


(1.7) 


Note that A{t) is orthogonal and hence satisfies A{t) ^ = A(t)*. We perform the change of unknown tp ijj 
defined by 

= ■ip{t,A{t)x). (1.8) 


This way, ip solves (1.2 1 if and only if ijj solves 


dti> = -Aijj - iV{t,A)ilj - 
where H is a time dependent potential given by 

V{t,±) = Vc{A{t )±), 

and the initial datum is 

• 0 ( 0 , x) = 0 o(x) = </ 5 o(x). 

For convenience, we shall denote explicitly the following time-dependent change of spatial variables: 


(1.9) 


Note that equation (1.91 satisfied by 0 is a standard nonlinear Schrodinger equation (113 with a cubic 
nonlinearity (k = 1 ) and a space and time dependent potential. 


1.2 General setting of the Cauchy problem 

The unboundedness of the potential function V makes the numerical analysis of the exponential methods 
difficult. Therefore, we modify the equation (113 by cutting off the potential V smoothly. Let us motivate 
this modification. In the context of Schrodinger equations with confining potentials, if the mass of the initial 
datum is essentially concentrated in a bounded set around the origin, then this mass localisation property is to 


be preserved by the evolution of (1.1), at least for reasonable times. Therefore, modifying the potential V out 


of a sufficiently large bounded set around the origin will not create huge errors in the solution, at least for not 
too long times. Let us introduce a smooth function y : K —> [0,1] such that 

Vx G [1 — (5/2, (5/2 — 1], x(x) = 1 and Vx G (— 00 , —(5/2) U ((5/2,-l-oo), x( 2 () = 0, 

where (5 2 is a given real number, chosen accordingly with the initial datum (po, the other physical parameters 

and the final computational time T > 0. We define the new potential function W for t G [0, T], x = (xi, • • • , XdY 
by letting 

d 

w( 0 x) = y( 0 x) J|x(xj). ( 1 . 10 ) 

i=i 

Although modifying the potential function as above changes deeply the physical situation, the fact that S is 
taken accordingly to ipQ and the physical parameters and for a finite time interval [0, T] gives some hope in the 
fact that the evolution of equation (1.1) with the potential function V and that of the same equation with V 


replaced with W starting from the same initial datum ipQ will be quite similar. 

To be more specific with the periodization of the problem, let us denote by the quotient K/((5Z). We 
consider the function w defined as (5-periodic in all directions such that for t G [0, T] and all x G satisfying 
\xj\ < (5/2, 1 < j < d, one has w(0x) = kF(t,x). The mapping t >->• w(t ,.) is smooth from [0,T] to as 

soon as cr > 0. In the following sections, we replace the continuous problem ([13 with its periodic counterpart 
as explained above, and we assume that (po = 0 o G for some a > d/2. 

We therefore get the following semilinear Cauchy problem in time: 


where T > 0, L = 


iA 


and 


i9t0(0x) - L'ip{t,x) = A‘u,(00)(x), (0x) G [0,T] X T/, 

0 (O,x) = 0 o(x), xGT/, 

A‘u)(00)(x) = -■jw(0x)0(t,x) - z/3|0p'^0(t,x). 


( 1 . 11 ) 


( 1 . 12 ) 
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Remark 1. The choice of the definition of the linear part L and the nonlinear part N from equation ([13 is 
somewhat arbitrary. One could also choose, for example 


i/s. 

L = — - iw{t,yi} 


and N{t,'tp){x.) = —ip\'tp\'^'^'tp{t,x.), 


but this would lead to a nonautonomous linear problem whose spectral properties are not as nice as that of the 
other case. However, in the case of a radially symmetric potential, one may get an autonomous linear part and 
use appropriate spectral methods on the whole space (see Eq. (2d) in M)- In this paper, we consider fairly 
general potentials without symmetry (see subsection |^.g| /or numerical experiments with ^ ly in (1-3)^ and 
we will indeed use the nice spectral properties of the operator L = iA/2 later. 

One can check that the operators of the form for a € K are isometric over The function 1 1 —> w{t, •) 

belongs to iL°’(T^)). Moreover, the nonlinear function iVy, satisfies a local Lipschitz condition: 

Lemma 2. For all T > 0, for all r > 0, there exists a constant C > 0 such that for all t G [0,T], ipi,ip 2 € 
iL'^(Tf) such that ||‘/5i||h^(t^) < r and \\t2\\h‘’{t)) < we have \\N„uit,ipi) - (^ 2 )||h^(t^) < C\\ipi - 

Remark 3. Note that this Lipschitz property is also true if one replaces N.^] with Nwit,-)^) = —iW(t, •)'!/'(■) — 
ifi\ijj\'^'^'tp{-) and the norms are taken in with a > d/2, but is no longer true with Ny = —iV{t, — 


We recover the fact that the Cauchy problem (1.11 1 is well-posed in H°'{Tg). 


1.3 Outline of the paper 

The outline of the paper is as follows: Section 2 is devoted to exponential Runge-Kutta methods. Firstly, we 
briefly describe the construction of these exponential methods, and recall what is the socalled underlying Runge- 
Kutta method [35], which is deflned using s collocation points. Secondly, we prove for the equation (1.11 1 that 
if the s collocation points are distinct, then the exponential Runge-Kutta method applied to Problem (1.11 1 
has order s. Moreover, we observe numerically that if we use Gauss collocation points, then we obtain order 2s 
(superconvergence). Section 3 is devoted to exponential integrators named Lawson methods. These methods are 
collocation methods on a new evolution equation, obtained from Problem (1.11) via another change of unknown. 


We show that Lawson methods applied to the simplified equation keep their classical order. In particular, the 
methods with s stages deflned with Gauss points have order 2s (superconvergence). Moreover, they preserve 
quadratic invariants up to round off errors. Section 4 is devoted to the comparison of our methods with other 
methods such as splitting methods as used in |5]. For our numerical experiments in dimension 1 and 2, the 
algorithm in time is supplemented with the fast Fourier transform (EFT) method in space, in a periodic domain. 
Finally we end the paper with conclusion and outlook. 


2 Exponential Runge-Kutta methods 


The main idea behind exponential integrators is to integrate exactly the linear part of the problem and then 
to use an appropriate approximation of the nonlinear part. Exponential Runge-Kutta (ERK) methods are 
particular exponential integrators. These methods have been derived and analysed for semi-linear parabolic 
Gauchy problems (see for example [35] for collocation methods, m for explicit methods and [35] for a survey 
on exponential integrators). They also have been used in |27] for solving linear and semi-linear Schrodinger 
Gauchy problems on the d-dimensional torus. They have been used to solve nonlinear Schrodinger equations 
for optical fibers (see the interaction picture method analyzed in m)- 


In this section, we introduce, analyse and use ERK methods to solve numerically equation (1.9l. 


2.1 Notations and description of the method 


In order to solve numerically problem ( |1.11 1 , we consider the following ERK methods of collocation type. We 
refer to |38] for a derivation of such methods for semi-linear problems based on variation-of-constants formula. 

Let T > 0 be the final computational time and {tn)o<n<M be a uniform subdivision of [0,r] with M +1 
points i.e. = nh with h = T/M. Let s S N* and Ci,--- ,Cs G [0,1] be given such that for all {i,j) G 
{1, • • • , s}^, Ci yf Cj if i j/ j. For some n G {0, • • • , M — 1}, we assume we have an approximation of the exact 
solution i({tn) of problem ( l.llj ) at time t„. Using an ERK method consists in computing an approximation 
V'n+l from ipn in the following way. First, we solve the nonlinear system of s equations 


V'n.fe = + h'^ak,t{hL)Ny,{tn + I < k < s, 


( 2 . 1 ) 


e=i 



















High order methods for nonlinear Schrodinger equations 


5 


where the unknowns are ,4’n,s and the coefficients {ak,i{hL)) 

operators on H°'{Tg) defined by 


1 ^ 
akAhL) = ^ 


with {Ci)i<e<s the Lagrange polynomials defined by 


c,(r)= n 




Cl — c,- 




are linear continuous 

( 2 . 2 ) 

(2.3) 


Then we compute ipn+i using 


'lpn-\-l = e^^ipn h^bk{hL)N^{tn -\- Ckh^'l/Jn^k), (2.4) 

fc=l 


where 


V/c e s}, bk{hL) = ^J Ck{a)da. 


(2.5) 


Remark 4. If we set L = 0 in formulae (2.21 and (2.5), then one recovers classical formulae defining Runge- 
Kutta collocation methods. The Runge-Kutta collocation method defined hy the corresponding coefficients is 
called the underlying Runge-Kutta method fMi- 


Note that in order to be able to compute fin+i from ifn, we have to precompute the coefficients a^A^A 
and bk{hL) for all fc, £ S {1, • • • , s}. We present an accurate and efficient way of precomputing these coefhcents 
in the next subsection. Of course, the spatial discretization of the operator L has to be specified. 


2.2 Precomputation of the coefficients 

It is well known that the Laplacian operator on iL'^(T^) is self-adjoint with eigenvalues 


UJp — 



{pI + ---+pA, 


P= (Pi,...,Pd) G 


We discretize the problem in space using a uniform grid with points and rely on FFT techniques. The 
eigenvalues of the discretized version of the periodic L = iAl2 operator are iujpj2 for 0 < |p|oo ^ K — 1, where 
IpIoo = rnax{|pi|,..., |pd|}. Since this operator is also self-adjoint, the computation of the discretized versions 
of the operators ak,i and bk amounts to the computation of the values of 

akA^hoJp/2)= -J^ ( 2 . 6 ) 

and ^ 

bk{ihujp/2) = ^ (2.7) 

We explain how one can compute these coefficients. These functions ak.i and b^ are holomorphic functions over 
C. After integration {e.g. using integration by parts) over (0, ft.) of the integrands defining these operators, we 
obtain holomorphic functions (say, of ftwp) with a removable singularity at the origin of the complex plane. For 
example, when the method has s = 2 stages, the coefficient ai i reads 


ai^i{ihu)p/2) = cf 


f _ jp^ft^^) -1-1- ih‘A^{c2 - Cl) 


(Cl - C2)((iCiftWp/2)2) 


( 2 . 8 ) 


Applying directly a formula such as (2.8) yields numerical instabilities in the computation of the ratio when 
ft|cLip| 1. Therefore, we use two different strategies for the computation of the coefficients (2.6) and (2.7). 


When ft-lwpl < 1/2, we use a discretized version of the Cauchy representation formula for a holomorphic function 

/ 


=A 


fA) 


do;, 


LJ — Z 


(2.9) 


(where C is the positively oriented unit circle and \z\ < 1), following the method presented in [l^l^T]. When 


ft|wp| > 1 / 2 , we simply evaluate directly formulas of the form ( 2 . 8 ). 
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Remark 5. For the evaluation of f{ihiOp) when h\u}p\ < 1/2, we use the trapezoidal quadrature rule to compute 
an approximation of the integral in (2.9). Since the unit circle is a smooth path in the complex plane, the 
function 1 1 —>■ ze®*/(e**)/(e** — ihujp) is a smooth 2'K-periodic function and hence the trapezoidal rule has infinite 
order. 

One can in fact compute these s x (s -|- 1) x coefficients independently and use parallel computing. 
Once these coefficients are computed, we can begin the time-stepping method and no additional computation 
is required provided the discretization parameters h and K are fixed. 


2.3 Result 


For the solution of the Cauchy problem (1.11), using an ERK method defined by the s points 0 < ci < 
Cs < 1, we have the following 


< 


Theorem 6. For all ifo G F[°'{Tg) and all T > 0 such that the exact solution of (1.11) is smooth over [0,T], 
there exists C, hg >0 such that for all h G (0,/ig), the ERK method defined by (2.1)-(2.4) starting from tjjQ is 
well-defined. Moreover, we have for all h G (0, ho) and n G N such that nh < T, 

W'f’iln) — ^ C'/l'* 

Proof. This is a consequence of Theorem 3.6 in m- Let us check the hypotheses of Theorem 3.6. Hypothesis 
3.1 is straightforward from (1.12). Hypothesis 3.2 follows from Lemma|^and the fact that is an algebra 

for a > d/2. Hypotheses 3.3 and 3.4 follow from our assumption on the temporal smoothness of the exact 
solution. □ 

Remark 7 . As mentioned in the introduction, the temporal constant C in Theorem is not uniform in the 
truncation parameter S. Note that our proof extends to general nonlinear Schrodinger equations with a smooth 
time-dependant potential on the whole space as soon as a > d/2 and the multiplication by V is a continuous 
operator from Fl'^ to itself. 

Remark 8. Even if we are not able to prove it, we observe in the numerical experiments that if the s collocation 
points of the ERK method are Gauss points, the method is of numerical order 2s (see subsection /.I). This 


indeed may be true when the methods are applied to the problem (1.11) on the torus T/ because of the periodic 


boundary conditions while it may not be true when the methods are applied to the problem (1.9) set on the whole 
space This is the case for example when one considers such methods applied to semilinear parabolic pro blem s 
This would indeed be another limitation of the reduction of the problem from to T/ (see subsection\l 


In any case, the proof of superconvergence for Gauss-ERK methods applied to problem (1.11) on the torus 


does not seem to be a straightforward adaptation of the classical result for ODEs. An option to prove it could be 
to use appropriate Taylor expansions of the exact solution of the problem, but the way one can order the terms 
appearing in the consistency error to obtain order 2s is really non trivial. 


In the following we call Gauss-ERK method an ERK method using Gauss collocation points. 


3 Lawson methods 


In [13 , Lawson considers the problem of designing some Runge-Kutta type methods for stiff ordinary differential 
equations. The idea is to perform a change of unknowns to transform the stiff system into a related nonstiff 
one. Then some basic Runge-Kutta method is applied to the related problem. The combination is termed a 
generalized Runge-Kutta method in m and often called Lawson method. 


The goal of this section is to describe the implementation of such methods on the problem (1.11) seen as an 


ordinary differential equation in time, to perform an analysis of the order of convergence of these methods as 
well as some of their preservation properties. 


3.1 Notations and description of the method 

Let us consider the equation (h 3. and set the following change of unknowns (so-called Lawson transformation), 

u{t,x.) = e~^*'i/{t,x.). (3.1) 


Then ip solves (1.11) if and only if u solves 


dtu{t,x.) = e e''"‘u(t,x)), (t, x) e [0, T] x Tf 

ii(0,x) = V'(0,x) = V'o(x), X e T^. 


(3.2) 
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Now one can apply a classical Runge-Kutta method to (3.21 seen as an ordinary differential equation in time. 
Assume {ak,i)i<k.i<s and {bk)i<k<s are a matrix and a vector of real entries. Set 


S 

VA: e s}, Ck ■■= 

e=i 


(3.3) 


and consider the s-stage classical Runge-Kutta method with Butcher tableau given by 


Cl 

fll.l 


Cs 

^S,l 

^S,S 


bi ■ 

• b. 


Assume that this Runge-Kutta method is of order at least 1. This means that the following condition is fulfilled: 


Y,bk = l. (3.4) 

k^l 

Applying this Runge-Kutta method to the problem ( |3.2[ ) defines a Lawson method: we compute an approxi¬ 
mation Un+i at time tn+i of the exact solution from an approximation at time tn by solving the system of 
s nonlinear equations (the unknowns being the {un,k)i<k<s)- 


un,k = Un + hY, + ah, , (3.5) 

e=i 

and then we compute Un+i through the formula 

S 

Un+1 =Un-\-h^ (^tn + Ckh, . (3.6) 

k^l 

Equivalently, the Lawson method for the unknowns 

V’n.z := and := (3.7) 

consists in solving the s nonlinear equations 

S 

4’n,k = e°'='"^V'n + (t„ -k ah, tpn.i), (3.8) 

ai 

and then computing 4’n+i through the formula 

S 

(t^ + Ckh.'ijjn.k) ■ (3.9) 

k^l 


We simply denote these relations ( |3.8| )-(3.9 1 by V'n+i = As for the ERK methods of (|^, the 

Runge-Kutta method defined by ak,i, bk is referred to as the underlying Runge-Kutta method. 

Note that, in view of Lemmathe Lawson method (3.7l-(3.9l is well defined in 77°’(T^) provided /i > 0 is 
small enough. 


3.2 Results 


In this section, we present some results on the Lawson method given by (3.8l-(3.9l. First of all, since equation 


(1.11) is time reversible, we give a sufficient condition for the Lawson method to be symmetric. We follow ideas 


developed in m. where the authors are solving an autonomous NLS equation, and we show that although 
equation (1.111 is non autonomous, the sufficient condition is quite similar. 


Theorem 9. Assume that the s-stage Runge-Kutta method defined by (ak,e)i<k,e<s o-nd (bk)i<k<s satisfies 


(3.41 so that it is of order at least 1. Assume that this method satisfies 


y(k,e)e sy 


O-s-kl-k.s+l-i + — 


(3.10) 


so that it is symmetric (see Theorem 2.3 in m)- Then the Lawson method defined by ( |3.8[ )-( [3i9| ) is also 
symmetric. 
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Proof. First of all, let us mention that the symmetry condition (3.101 gives 

be = bs+i-e, 1 < i < s. 


Moreover, summing (3.101 over i and using (3.41 we have 


1 Cf^ — Cs-|_i_fc, 1 ^ ^ S. 


(3.11) 


(3.12) 


The adjoint of the method is by definition The relation ifn+i = is 


equivalent to ipn = This corresponds to exchanging tn with and h with —h in (3.8)-( [3^ . 

This leads to 


^n,k = e (^tn+l — Ceh,'tpn,ej 

i=l 

S 




Extracting ifn+i from (3.14) gives 

S 

Ipn+1 = e'‘^'0n + ^ X! (tn+1 - Ckh, ■ 




Plugging this expression into (3.131, we get 


s 

tpn.k = + h'^{be - (tn+i - ceh,-tl>n,e'^ ■ 


Using (3.10), we infer 

S 

4’n,k = -(- (t„+i - ceh,i>n,ej 


e=i 


Setting ipn.k = 'ifn.s+i-k and reordering the sums in (3.151-(3.161, we obtain 

S 

^u,k = + h ^ (t„ + (1 - c,+i_,)h, , 


(3.13) 

(3.14) 

(3.15) 


(3.16) 


■ifn+l = e^V^n + ^ X! (tn + (1 “ Cs+l-fc)h, ■0„,fc) . 


k^l 


Using (3.111-( 3.121 we conclude that 

s 

'^n,k = + h'^ (t„ -h ceh, -ipn.e) 


1pn+l = e^^-lpn + h^ bkt''^ Vu, (tn + Ckh, ^n,k) ■ 


fc=l 


This proves that the Lawson method is symmetric. 


□ 


If the underlying Runge-Kutta method preserves quadratic invariants in the sense that it satisfies the Cooper 
condition, then so does the associated Lawson method. 


Theorem 10. Assume that the underlying Runge-Kutta method satisfies (3.4) so that it is of order at least 1. 
Assume that it satisfies the Cooper condition, 


bko-k.t + beae^k = bkbe, V 1 < fc,^ < s. 


(3.17) 


so that it preserves quadratic invariants. Then the Lawson method (3.8)-(3.9) preserves the L‘^-norm: 

ll^n||L 2 (T^) = II'0 o||l2(t^)i V n > 0. 
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Proof. Since the evolution equation (1.111 preserves the L^-norm, and so does the change of variables (3.1), the 
evolution equation (3.21 also preserves the L^-norm. Therefore t i—>■ x) pdx is a quadratic invariant of 


equation (3.21. The Cooper condition guarantees that the Runge-Kutta method (3.5l-(3.6l preserves quadratic 
invariants. Hence for all n, J'.j.d |u„(x)pdx = |uo(x)pdx and then |■^/)„(x)pdx = j^|V'o(x)pdx. □ 

such that the under- 


Definition 11. We eall Gauss-Lawson method any Lawson method of the form 
lying Runge-Kutta method is a Gauss eolloeation method (see section II.1.3 in m)- 


Corollary 12. A Gauss-Lawson method applied to (|1.11| is symmetric and preserves the L^-norm. 


Proof. The underlying Runge-Kutta method is a collocation method at Gauss points and therefore it is sym¬ 


metric (see Corollary 2.2, Chapter V in [M|)- Hence the Gauss-Lawson method applied to (1.111 is symmetric 
by Theorem]^ Moreover the underlying Runge-Kutta method satisfies the Cooper condition (3.17) (see exercise 
5 in chapter IV in [M])- Therefore the Gauss-Lawson method preserves the L^-norm by (10). □ 


We now want to prove that the Gauss-Lawson method with s stages applied to (1.11) has order 2s in 11°'(Tg) 
for a > dl2 (see Theorem 14). Our strategy is the following. First, we consider an equivalent autonomous form 
of (1.11). Second,we show that applying a Lawson method to this autonomous form is essentially the same as 


applying the method to (1.11) directly. Third, we rely on an Alekseev-Grobner lemma for autonomous systems 


which provides a representation of the error that allows us to conclude that, since Gauss-Lawson methods are 
collocation methods, they have order 2s. 

Let us set 

^ ^ TB NX 

.u{t) 


Uit) = 


X H'^(T^), 


so that (1.11) reads 


with 


We have the following useful 


F 


-U{t) = F{U{t)), 




(3.18) 

(3.19) 


Lemma 13. Let us fix uq G II°(Tg). For all h > 0 sufficiently small, we apply the Lawson method (3.5)-(3.6) 
and denote by Un the corresponding numerical values. Similarly, we start with Uq = ( 0 , mo )*> we apply the same 
method to the Gauchy problem ( |3.18[ ) with intial datum (7(0) = Uq and we denote by [/„ the corresponding 
numerical values. We have for all n G N such that nh < T, 


U„ = 


nh 

u. 


Proof. Since the function Ny, satisfies a local Lipschitz condition (see Lemmaj^, so does F (see (3.19)). Hence, 
the Lawson methods are well defined locally for h > 0 small enough. We perform the proof by induction. The 
relation is true for n = 0. Assume it holds for some n S N such that {n l)h < T. The definition of the 


coefficients {ck)i<k<s (see (3.3)) and the first component of the function F ensure that for all fc G {1,..., s}, the 
first coefficient of Un^k is t^ k = + Ckh. Therefore, in view of (3.5), the {Un^k)i<k<s = {n-h + Ckh,Un^kY are 

the unique local solutions for the Lawson inner problem in U. Similarly, the consistency relation (|3.4[) ensures 


that the first component of (7„+i is (n -|- l)h. Hence, in view of (3.6), we have that the second component of 
Un+i is actually Un+i. 

We are now able to state 


□ 


Theorem 14. Assume uq G II°(Tg) is given and T > 0 is chosen sueh that the exact solution 1 1 —>■ u{t) of the 


Gauehy problem (3.2) is defined over [0,T]. There exist constants C > 0 and Hq > 0 such that the eorresponding 
numerical approximations provided by the Gauss-Lawson method with s stages Un satisfy 

Vft. G (0, ho), Vn G N s.t. 0 < nh < T, IWi'tn) — Un\\H<^(T‘^) ^ Ch^'^. 


Remark 15. This implies that the exact solution if of (1.11) also satisfies 


Wh G (0, ho), Vn G N s.t. 0 < nh < T, Wf’i'tn) — f’nW h<^ (T'^) ^ Ch^^, 


where ipn cire the numerical approximations of ip by the Gauss-Lawson method (3.8)-(3.9). 
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Remark 16. As mentioned in the introduction, the temporal constant C in Theorem 14 is not uniform in the 
truncation parameter S. 


Proof. Thanks to Lemma it is sufficient to prove that the Gauss-Lawson method applied to the Cauchy 
problem (3.181 with initial datum Uq = (0,uo) is of order 2s. We chose hg > 0 sufficiently small to ensure 
that the method is well-defined for all h € (0, ho). We follow the proof of [33] and [H] for collocation methods 
applied to ODEs. Indeed, with our hypotheses (cr > d/2,K € N) the vector field F defined in ( 3.19[ ) is 
C°“(K X lI'^(Tg),R X iL'^(T^)) and we are dealing with an autonomous ODE in the Banach space K x 
with smooth vector field. Hence, solutions of ( 3.18[ ) are smooth functions with values in K x Therefore, 

the variational evolution equation for Y = 

Y'{t) = F'{U{t))Y{t), 


with V(0) = (l,IdH<^(Td)) has a right-hand side linear operator t ^ F'{U{t)) S C°°(/,£(K x iL‘^(T^))), where 
I is an interval where U is defined and £(K x iL'^(T^)) stands for the space of linear continuous operators of 
M X iL'^(T^). First, we recall that the mapping C/„ '—>■ Un+i can be defined through a collocation problem. 
The same problem allows us to express the consistency error: find P a polynomial of degree at most s (with 
coefficients in K x such that 


r P{tn) = Uitn), 

\ P'{t„ + Ckh) = F{P{tn + Ckh)), for I < k < s, 

and set Un+i = P{tn+i), so that the consistency error reads Un+i — U{tn+i). Second, define the defect d as 

d{t,P{t)) = P'{t)-F{P{t)). 


Then, we use the Alekseev-Grobner lemma (Theorem 2 of |13]), which ensures that there exists a smooth 
function g (corresponding to 8280 in US) from [0,T] x (M x iL'^(T^)) to the linear continuous operators of 
R X iL°’(T^) such that the following identity holds in iL°’(T^): for all t € [0,/i], 

P{tn + t)-U{tn + t) = J g{tn + t-T, P{T))d{T, P{T))dT. 


Setting t = h in this formula and using the fact that g is smooth and all the derivatives of t 1 —>■ d{t, P{t)) are 
bounded on [0,/i] independently of h (see Lemma 1.6 of Ghapter 2 in |M|), we infer that the consistency error 
in is of order Hence the global error is of order hf^ via a classical discrete Gronwall lemma. □ 


4 Numerical experiments 

In this section, we make some numerical experiments and show the numerical efficiency of the methods. First, 
we consider one dimensional problems and compare some ERK methods and Gauss-Lawson methods with the 
classical splitting methods. Second, we compare the same kinds of methods applied to an actual 2D problem 
with a rotating BEG. As expected in the theoretical results above , the numerical behavior does not actually 
depend much on the dimension as is illustrated in subsections |4.1| and |4.3| 


4.1 One dimensional example: the one-dimensional cubic NLS equation 

We provide in this subsection some numerical experiments to show the efficiency of the Gauss-ERK and the 
Gauss-Lawson methods compared with the traditional splitting methods. We use here splitting methods re¬ 
spectively of order I, 2, 4 and 6 [Tl]. In order to present the usual splitting schemes, we have to define the 
operators Sl and associated respectively to the evolution of the equations 


dtu{t,x) = Lu{t,x) dtv{t,x) = Vu,w(t, x). 


where L = iA/2 and N^, is defined in (I.I2). The operators satisfy by definition the following relations involving 
the exact solutions of the associated equations: 


u{t + h,x) = SL{h)u{t,x) and v(t + h,x) = SN„ih)v{t,x). 

As explained for example in a splitting idea for building a splitting scheme of even order consists in 
approximating the continous fiow associated to dtC{t,x) = Lf{t,x) + Ni„{t, (4{t,x)) by a composition of operators 
Sl and 5'^?^ of the form 

Sn„ {aih)SL{bih) ■ ■ ■ (orh)SLiKh)S(or+ih), 
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Order 2 

Order 4 

Order 6 

( 

II 

II 

oi = 0.0502627644003922 

6 i = l 

02 = 1/2-9 

02 = 0.413514300428344 


bi = 29 

03 = 0.0450798897943977 


62 = 1 - 40 

04 = -0.188054853819569 


9 = (2-h2i/3 + 2 -i/3)/6 

05 = 0.541960678450780 

Og = 1 — 2 (oi -|- 02 -l- 03 -|- 04 -|- 05 ) 

61 = 0.148816447901042 

62 = -0.132385865767784 

6 3 = 0.067307604692185 

64 = 0.432666402578175 

65 = 0.5 — (5i -l- 62 + ^3 + bi) 


Table 1: Coefficient of splitting methods. 


with r > 1 and for all ar+ 2 -e = de and br+i-e = bi. The coefficients of the splitting methods we implemented 
are given in Table ^ Let us mention that alternative ways of deriving splitting methods for such problems 
exist, some of which authorize the use of complex coefficients. For example, one can use splitting methods such 
as (5.1) in [TT] , where only one set of coefficients has nonzero imaginary part and the other set can be used to 
solve the linear part of a Schrodinger equation (see also [I^). 

We compute the numerical solution to the one-dimensional cubic NLS equation 

dtip = idlip + (t,x) G [0,r] X M. (4.1) 


An exact solution for (t, x) G [0, T] x K is given by the soliton formula 


‘4^ex (tj x'j 


—sech((a; — xq) — ct) exp 



xo) — ct 
~2 



(4.2) 


For the numerical approximation of this solution, as explained in subsection |1.2| we take a periodic finite interval 
{xi,Xr) of big enough length, and we discretize the space operators using Fourier spectral approximation. We 
choose the spatial mesh size k = Ax > 0 with k = {xr — xg)/M with M = 2^, P G N*. The time step is denoted 
by = At and h = T/Nt for some Nt G N*. The grid points and the discrete times are 


Xj:=Xi+jk, tn:=nh, j = 0,l,---,M, n = 0,1,- 




T- 


Let ip'j be the approximation of ip{tn,Xj). Since we discretize (4.11 by the Fourier spectral method, and its 
Fourier transform satisfy the following relations: 


M/2-1 

m=-M/2 




J = 0, 


,M-1, 


and 


M-l 




where fj,m = for all m = — • • • , ^ — 1. Let us define the discrete gradient operator 


3=0 

2 ’ ■ ■ ■ I 2 

= 'i-k-mVm, V G 

Let us denote by H^ the projection operator 

Hfc : C^{[x^,Xr\,C) 




{‘piXj)) 


>0<3<M-1 


We define the discrete P' norm on as 


M-l \ 

, uGC^, r>l. 

3=0 ) 

Using these definitions, we consider the following errors: 

^P.h — sup ||n/j(l/:ea;(tnj *)) (l^7)i|L2 ? 

nG{ 0 ,---.Af} 



(4.3) 
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-M,h — 


sup 

,N} 


||n/j; {ipex (j'Tl 


\(Q. 


-WmjL. /l|nfc(V'ex(o, 


\p 


We also define the discrete energy: 

Ek{v)=\\\^kv\\%-l\\v\\%, 

the energy conservation is seen through the following relative error 

£E,h= sup \Ek{Tik{i’ex{tn,-))) - 
,N} 


(4.4) 


(4.5) 


For all the following simulations, we consider the computational domain (xe,Xr) = (—15,15) and the final 
time T = 5. The experiments are performed with the discretization parameters P = 10 and various time steps h 
and the chosen physical parameters are q = 8, a = q^/16 and c = 0.5. We provide for all methods the evolution 
of £p.h, ^M.h and £E,h for various time steps h. All of our methods being implicit, it is necessary to solve a 
nonlinear problem at each time step for both ERK and Lawson method. This is performed here through a fixed 
point algorithm. It is therefore important to show the efficiency of the schemes to plot the evolution of the CPU 
time with respect to the error £p^h- 


The legends of the figures are respectively given in Figures 1(a) and |l(b)] for ERK and Lawson methods 
(we recall that s denotes the number of collocation points) and splitting schemes. We showed in the previous 
sections that if one uses arbitrary distinct collocation points, the order of ERK and Lawson methods are h®. 
Therefore, in all figures, curves associated to s = 1, 2 and 4 respectively correspond to splitting of order 1, 2 
and 4. If we would like to compare the splitting of order 6 to ERK and Lawson methods, we would have to 
consider s = 6 collocation points. 


s = 1 


-0- s = 2 


—^— Splitting 0(1) 

s = 3 


—e— Splitting 0(2) 

—s = 4 


—0— Splitting 0(4) 

s = 5 


□ Splitting 0(6) 


(a) ERK and (b) Splitting schemes 

Lawson meth¬ 
ods 


Figure 1: Legends 


We gather the evolution of £p^h for various schemes in Figure]^ Note that since we use Gauss collocation 
points, we clearly obtain numerically an order 2s for both Gauss-ERK and Gauss-Lawson methods with s 
stages. This is predicted for Gauss-Lawson schemes thanks to Theorem Let us remark that there is a 
saturation phenomenon for high s: the phase error £p^h ceases to decrease when h decreases. This is due to 
the fact that the exact solution defined in (4.2) do es no t satisfy exactly the boundary conditions (see the 
transformation of the problem from to in section 1.2). This can be reduced either by increasing the size 


5 of the computational domain or by changing the periodic boundary conditions to, for example, homogeneous 
Dirichlet ones (in this case, one can replace EFT with the discrete sine transform). If we compare Figures [2 (a) | 
and |2(b)] it is noticeable that the constants of error are better for Gauss-ERK methods. 

The preservations of mass and energy are fundamental when one deals with dispersive equations. Let us 
first discuss mass preservation. The figure presented in Figure]^ shows the evolution of £M,h with respect to the 
time step for several Gauss-ERK methods. Indeed, Gauss-Lawson methods preserve mass up to round-off (see 
Theorem 10) and so do splitting methods. It is noticeable, however, that, on the numerical example (4.1)-(4.2) 
a Gauss-ERK method with at least two stages achieves round-off in mass preservation for time steps below 
10-2. 


Goncerning the energy conservation, none of the presented methods is able to handle it properly. However, 
it is clear in Figure]^ which displays the evolution of £E,h as a function of the time step for all schemes, that, 
with respect to energy preservation, Gauss-ERK methods are of better quality than Gauss-Lawson ones which 
are themselves better than splitting schemes. 

Finally, Figure [^shows the evolution of the GPU time with respect to the £p^h error. Theses figures are very 
interesting since one could think that implicit methods are more costly compared to splitting schemes. For a 
given error £p^h, the GPU time is clearly lower for Gauss-ERK schemes. For a “big” phase error £p^h > IQ-^, 
the splitting schemes are less costly than Gauss-Lawson ones. The situation is reversed when one is interested 
in “small” phase errors £p^h < IQ-^ and Gauss-Lawson schemes have to be recommended. For both ERK and 
Lawsons schemes, it is not necessarily interesting to use a high number of collocation points. 
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(a) Gauss-ERK (b) Gauss-Lawson 



(c) Splitting 


Figure 2: Evolution of £p^h error 


with respect to the time step for problem (4.11 (see Figure [^for legends) 



Figure 3: Evolution of error with respect to the time step for problem (4.11 using Gauss-ERK methods 
(see Figure for legends) 


Another interesting question is the behaviour of conserved quantities (mass, energy) by the methods over 
long times. For this numerical study, the final time is set to 5000 and to avoid any hazardous behaviours due 
to boundary conditions, we choose as initial condition ip{0,x) = |sin(a;)| and x G [—7r,7r). With this choice, 
the datum belongs to The number of Fourier modes is and we take h = 10“^. We are interested 

in the quantities £M,h{tn) defined in (4.4| and £E,h{tn) defined in ( |4.5| ). We plot in Figure [^the evolution of 
£M,h{tn) for both Gauss-ERK and Gauss-Lawson methods. The L^-norm conservation holds up to roundoff for 
Gauss-Lawson methods as predicted by Theorem 10 Almost preservation of the L^-norm over long times seems 
to hold for Gauss-ERK methods. In contrast, the evolution of £E,h{tn) (Figure]^ shows that the Gauss-Lawson 
scheme breaks down in such a situation. The Gauss-ERK scheme seems to almost preserve the energy and 
is clearly superior, with respect to energy preservation, to the Lawson scheme for long time simulations. The 
break down phenomenon of the Lawson scheme can be reduced by selecting a smaller time step {h = 10“^, see 
Figure]^. However, the numerical cost is dramatically more important compared to a Gauss-ERK scheme. 

In conclusion, even if the Gauss-ERK schemes do not preserve theoretically neither the norm nor the 
energy, they allow us to preserve them numerically for reasonably small time step h even in long time simulations 
and clearly are the best methods for one dimensional simulations of the cubic NLS equations. 


4.2 One dimensional example: a NLS equation with time-dependent potential 
and cubic-plus-quintic nonlinearity 
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(a) Gauss-ERK (b) Gauss-Lawson 



(c) Splitting 


Figure 4: Evolution of £E,h error with respect to the time step for problem (4.11 (see Figure [^for legends) 




(a) Gauss-ERK (b) Gauss-Lawson 



£p,h 

(c) Splitting 


Figure 5: Evolution of the CPU time with respect to £p^h error 


for problem (4.11 (see Figure]^ for legends) 


We consider here the following one-dimensional NLSE (see e.g. [TU]) 

idttp = -dl'ip + Vtp + Gi\tp\‘^4’ + G2\^fip, (4.6) 

with a cubic-plus-quintic nonlinearity and an external time-dependent potential given by V{t, x) = |tLi^ cos(a;t). 
The splitting schemes have to be applied carefully for this equation since the step associated to the potential 
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Figure 6: Evolution of SM,h{tn) (figures (a) and (b)) and SE,h{in) (figures (c) and (d)), /i = 10 ^ 


10 ® 


10 ° 


Ki 10- 


10 - 


10 - 1 ® 


-Lawson 5 = 1 

-Lawson s = 2 

-Lawson s = 3 



0 


2500 5000 

t 


Figure 7: Evolution of SE,h{tn), h = 10 Gauss-Lawson 


term leads to a non-autonomous problem. However, the Gauss-ERK and Gauss-Lawson methods can be applied 
straightforwardly for this test case. 

We know in this case an exact solution (the bright soliton) given by 


exp 
= V - 



- ‘^x sin{ujt + Po) - sin(2a;< -I- 2/3o) - Ect^ ^ 

(Vl — bcosh{2^/—Ec{x — cos(a;f))) -I- l)^^^ 
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The final time T is chosen to be 5. The parameters of the reference solution are fixed to 

Gi = —2, G 2 = -, w = 2, Ec = —1, /3o = 0, T] = ^ = —16 . 

The spatial computational domain is T) :=] — 32, 32[, with = 2^^ Fourier modes. For this example, the energy 
is not conserved. We therefore present in Figure]^ only the evolution of £p^h and £M,h errors. We only draw 
the error curves for s < 4 for Gauss-ERK and Gauss-Lawson methods since the curves for s = 5 are identical 


to the ones for s = 4. and since Gauss-Lawson methods preserve the mass up to round-off (see Theorem 101. 


Like for the cubic NLS equation in subsection 4.1 we have good decay rate of the £p^h curves. 




(a) Gauss-ERK (b) Gauss-Lawson 




(c) Gauss-ERK (d) Gauss-Lawson 


Figure 8: Evolution of £p^h (figures (a) and (b)) and 
step for problem (4.6 1 (see Figure for legends) 


£M,h (figures (c) and (d)) errors with respect to the time 


4.3 Two dimensional simulations 

The first 2D example concerns the equation 


StV'= (i,x) e [o,T] X 


(4.7) 


which is a nonlinear Schrodinger equation without confinement or rotation. We look for a solution of the form 

V'(t,x) = e*‘0(x). 


where 0 solves 


A0-k|0p0 = 0, 


with lim||x||_>.oo 0(x) = 0. Since we do not have access to an exact solution 0 to this problem, we generate it 
through a classical shooting point method |26]. Since the decay of this kind of solutions is slow, we consider 
(x, y) £ [—38, 38]^ with P = 9. 

As for one dimensional simulations, we evaluate several errors for the different methods. We generalize the 


error functions £pm, £M,h and £E,h defined in (4.3l, (4.4) and (4.5) to two dimensional cases. We only present 
here simulations for ERK and Lawson methods with s = 1, 2 and 3 since we have seen on one dimensional 
experiments that we do not gain much more precision for s > 3. We also restrict ourselves to the presentation of 
the £p^h error displayed in Figure |^since the £M,h and £E,h errors lead to similar results to the ones presented 
in Figures]^ and 
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(a) Gauss-ERK (b) Gauss-Lawson 



(c) Splitting 


Figure 9: Evolution of Ep^h error with respect to the time step for problem (4.71 (see Figure for legends 


We also propose in this subsection to illustrate the efficiency and versatility of the methods on the simulation 
of a soliton for the two-dimensional Schrodinger equation and the simulation of the evolution of a rotating Bose- 
Einstein condensate. Thereby, our second 2D numerical example corresponds to a rotating BEC modelled by 


equation (1.2). We reproduce here with the Gauss-ERK method of order 6 the simulations realized in [5] for 
(1.2), see Figure 10 The parameters are (3 = 1000, El = 0.9, and the potential is (1.3) with = 1.05 and 
7 j, = 0.95. The computational domain is (—16,16)^ with 2® Fourier modes in each direction. The time step 
is h = 10“^ with a final time T = 7. The initial datum is the ground state of the stationary equation and 
was generated using Matlab toolbox GPELatQ We recover the same behaviour and we get good conservation 
properties. 


Conclusion 

We presented ERK and Lawson methods that allow us to compute numerical solutions of NLS equations 
modelling a rotating BEG in a fairly general setting. This procedure allows us to derive neatly high order 
methods (in contrast to finding coefficients for high order splitting methods with small error constants, for 
example). We proved finite-time convergence in Sobolev norms for these methods in a simplified framework. 
We compared the numerical results provided by these methods to that obtained via classical splitting methods 
in several configurations: ID problems, 2D problems without confinement or rotation, 2D realistic problems 
with rotation and confinement. It is noteworthy that all the methods presented in this paper allow to deal with 
non-autonomous problems (no matter whether the lack of autonomy comes from the physical situation or a 
change of unknown). 

When it comes to finite time accuracy, Gauss-ERK methods outperformed Lawson methods and splitting 
methods in all cases, since they have very low error constants numerically. Moreover, even if they do not preserve 
mass exactly (up to roundoff errors), their relative error on the mass is comparable to that of Gauss-Lawson and 
splitting methods (which preserve the mass up to roundoff errors) for reasonably small time steps (see Figure]^ 
for time steps of order 10“^). This can be explained by the accumulation of roundoff errors when h gets smaller 
for the methods preserving mass exactly. 

When it comes to computational times, Gauss-ERK methods, though implicit, also outperformed Lawson 
and splitting methods (see Figure]^ on the ID example. To be completely fair, one must say that ERK methods 
(and, to a lesser extent, Lawson methods) require the precomputation of some coefficients. This computation 
only needs to be carried out once, before one starts the time stepping procedure, and can be parallelized if 
necessary. 

^http://gpelab.math.cnrs.fr 
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t=0 t=1.5 





Figure 10: Contour plots of the density function |(/?(t,x)p in a rotating BEC. 


With respect to average or long time behaviour, in contrast to splitting methods, ERK and Lawson methods 
do not show resonances phenomena UHl El- 
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